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A  computational  study  of  the  electrochemical  hydrodynamic  process  in  an  alkaline  fuel  cell  was  con¬ 
ducted.  The  computation  relaxed  the  ideal  solution  assumption,  accounted  for  thermodynamic  solubility 
of  the  reactants,  and  allowed  for  property  variations  due  to  temperature  and  concentration  effects.  The 
results  showed  that  the  ideal  solution  assumption  is  not  adequate  for  calculation  of  the  transport  process 
of  the  concentrated  electrolyte  considered,  7  M.  The  ideal  solution  formulation  resulted  in  a  lower  lim¬ 
iting  current  density  condition  by  about  50%  than  that  predicted  by  the  non-ideal  solution  formulation. 
The  study  also  showed  that  the  thermal  condition  is  important  to  the  calculation  of  the  limiting  current 
density  condition.  The  calculated  limiting  current  density  increased  by  about  30%  when  the  boundary 
condition  was  changed  from  isothermal  to  adiabatic.  The  computational  results  suggest  that  maintain¬ 
ing  a  uniform  KOH  concentration  in  the  electrolyte  (for  example,  at  design  point  of  7  M)  be  an  effective 
measure  to  increase  the  limiting  current  density  condition. 

©  201 1  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Stationary  fuel  cell  is  a  viable  option  for  electrical  power  gener¬ 
ation  from  renewable  sources.  For  example,  H2  produced  by  solar 
or  wind  turbine  based  electrolysis  can  be  used  to  fuel  the  fuel  cells. 
Among  the  fuel  cell  platforms,  alkaline  fuel  cell  (AFC)  remains  a 
good  prospect  as  it  has  high  reliability  and  it  can  use  non-precious 
metal  for  the  electrodes  [1-6,23].  However,  to  be  competitive  with 
other  power  sources,  the  power  density  (as  well  as  the  cost  of 
manufacturing,  operation  and  maintenance)  needs  to  be  improved. 
A  CFD  based  study  of  the  electrochemical  hydrodynamic  process 
that  sets  the  limiting  current  density  condition  is  conducted.  The 
purpose  of  this  study  is  to  gain  a  better  understanding  of  the  mech¬ 
anisms  critical  to  setting  the  limiting  current  density  condition, 
and  to  demonstrate  the  application  of  CFD  based  simulation  to  aid 
design  consideration. 

The  CFD  formation  accounts  for  the  electrochemical  reaction, 
charge  and  species  transport,  and  thermodynamics  of  gas  solubil¬ 
ity  in  the  electrolyte.  The  mathematical  models  that  describe  the 
transport  processes  are  available  in  the  literature,  e.g.,  see  [7-13]. 
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However,  most  transport  models  were  based  upon  the  assumption 
of  ideal  solution  for  the  calculation  of  the  transport  properties.  Error 
associated  with  the  ideal  solution  assumption  was  not  assessed. 
Furthermore,  most  CFD  studies  reported  in  literature  assumed  an 
isothermal  system.  The  temperature  effects  were  not  evaluated.  In 
this  paper,  the  ideal  solution  assumption  was  relaxed  and  the  tem¬ 
perature  effects  were  examined.  The  objectives  of  this  paper  are 
(a)  to  examine  the  effects  on  CFD  calculation  due  to  the  ideal  solu¬ 
tion  and  isothermal  assumptions,  and  (b)  to  gain  insights  into  the 
mechanisms  that  are  critical  to  setting  the  limiting  current  density 
condition. 

2.  Formulation 

2.  \ .  Transport  equations 

Major  assumptions  invoked  are  (a)  Newtonian  fluid,  (b)  negli¬ 
gible  pressure  effects  on  enthalpy,  (c)  negligible  magnetic  effects 
due  to  electrical  field,  (d)  negligible  pressure  and  temperature 
effects  on  diffusion  transport,  (e)  negligible  viscous  dissipation, 
pressure  work,  and  Dufour  effects  in  energy  equation,  (f)  no  homo¬ 
geneous  chemical  reactions,  (g)  electroneutrality,  (h)  homogeneous 
and  continuous  media  for  gas  and  liquid  phase,  (i)  local  thermal 
equilibrium,  (j)  vapor-liquid  equilibrium  at  the  gas-liquid  inter¬ 
face,  and  (k)  negligible  electrical  resistance  of  the  electrodes.  The 
governing  equations  describing  the  transport  of  mass,  momentum, 
species,  energy  and  charge  balance  are  summarized  in  Table  1.  The 


4924 


G.  Zhou  et  al.  /  Journal  of  Power  Sources  196(2011)  4923-4933 


Nomenclature 

ak  activity  of  species  k 

dpy  specific  interfacial  area  between  phase  j3  (liquid) 

and  phase  y  (gas)  (m-1 ) 

ap0  specific  interfacial  area  between  phase  (liquid) 
and  phase  o  (solid)  (m-1 ) 

Ck  molar  concentration  of  species  k  (kmol  m-3 ) 

Cp  specific  heat  (J  kg-1 I<-1 ) 

Dk  diffusion  coefficient  of  species  k  into  the  solu¬ 
tion/mixture  (m2  s-1 ) 

Dkj  multicomponent  Maxwell-Stefan  binary  diffusion 

coefficient  (m2  s-1) 
df  fiber  diameter  (m) 

E  electric  potential  (V) 

F  Faraday  constant  (C  kmol-1 ) 

g  gravitational  acceleration  (ms-2) 

H  Henry’s  constant  (kmol  m-3  atm-1 ) 

i  current  density  (A  m-2 ) 

z'o  exchange  current  density  (A  m-2 ) 

i  vector  quantity  of  current  density  (A  m-2 ) 

in  current  density  at  solid-liquid  interface  (A  m-2 ) 

jk  vector  quantity  of  mass  flux  of  species  k  (kg  m-2  s-1 ) 

I<  absolute  permeability  (m2) 

kr  relative  permeability 

kK  Kozeny  constant 

kSCx  Setschenow  salt  effect  parameter 

Mk  molecular  weight  of  species  k  (kg  kmol-1 ) 

m  molality  (mol  kg-1 ) 

n  number  of  electron  transfer 

p  pressure  (Pa) 

pc  capillary  pressure  (Pa) 

ps  saturation  vapor  pressure  (Pa) 

qk  stoichiometric  coefficient  of  species  k 

Ru  universal  gas  constant  (J  kmol-1  K-1 ) 

5  entropy  (J  kmol-1  K-1 ) 

Sr  reduced  phase  saturation 

s  saturation 

Sim  immobile  saturation 

sk  stoichiometric  coefficient  of  species  k 

T  temperature  (K) 

t  time  (s) 

tk  transference  number  of  species  k 

U  open  cell  potential  (V) 

uk  mobility  of  species  k  (m  s-1  kmol  N-1 ) 
v  velocity  vector  (ms-1) 

xk  mole  fraction  of  species  k 

Yk  mass  fraction  of  species  k 

zk  charge  number  of  species  k 

Greek 

a  transfer  coefficient 

8  film  thickness  (m) 

y±tc  mean  activity  coefficient  based  on  molarity 

s  porosity 

p  local  overpotential  (V) 

k  electrical  conductivity  (S  m-1 ) 

kd  diffusion  conductivity  (A  m-1 ) 

A  thermal  conductivity  (W  m-1  K-1 ) 

A°  limiting  ionic  equivalent  conductance  of  ion  k 
(S  m2  kmol-1 ) 

fi  dynamic  viscosity  (kg  s-1  m-1 ) 

pLk  chemical  potential  of  species  k  (J  kmol-1 ) 

v  stoichiometric  coefficient 


vk  number  of  cation  or  anion  produced  by  the  dissoci¬ 
ating  electrolyte 
0C  contact  angle 

p  density  (kg  m-3) 

a  surface  tension  (Nm-1) 

r  normal  or  shear  stress  (Nm-2) 

0  electric  potential  at  electrolyte  phase  (V) 

Subscripts 
a  anode 

c  cathode 

e  electrolyte 

8  gas 

ij,k  species 

l  liquid 

m  mixture 

nw  non-wetting 

T  total 

w  wetting 

P  liquid  phase 

y  gas  phase 

o  solid  phase 

0  standard  state,  pressure  at  1  atm  or  solvent  or  refer¬ 

ence  value 
+  cation 

-  anion 

Superscripts 
a  anode 

c  cathode 

eff  effective 

h  energy 

i  charge 

m  mass 

r  reference  state 

j5  liquid  phase 

y  gas  phase 

o  solid  phase 

0  standard  state,  pressure  at  1  atm,  or  sol¬ 

vent/reference  state 


respective  source  terms  are  given  in  Table  2.  The  source  terms 
account  for  (a)  mass  and  species  addition  or  removal  due  to  gas 
solubility  and  H20  phase  change,  (b)  momentum  exchange  due  to 
Darcian  flow  in  porous  electrodes,  (c)  energy  source  or  sink  terms 
due  to  heat  and  entropy  generation  at  catalyst  layer,  Joule  heat¬ 
ing  in  separator  and  catalyst  layer,  and  latent  heat  of  H20  phase 
change,  and  (d)  charge  generation  at  catalyst  layer.  The  constitutive 
equations  describing  the  reaction  rate  (Butler-Volmer  equation), 
capillary  pressure  (Leverett-J  function),  and  physical  properties  are 
given  in  Table  3.  Conservation  equation  of  the  liquid-phase  species 
is  solved  only  for  OH-.  The  mass  fraction  of  K+  is  calculated  from 
electroneutrality. 

2.2.  Electrochemical  reaction 

One-step  electrochemical  reactions  are  assumed  for  anode  and 
cathode  electrodes, 

Anode: 


H2  +20H-  ->  2H20  +  2e~ 


(1) 
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I 

KOH 


1:  Anode  Gas  Channel 
2:  Anode  Gas  Diffusin  Layer 
3:  Anode  Catalyst  Layer 
4:  Separator  Layer 
5:  Cathode  Catalyst  Layer 
6:  Cathode  Gas  Diffusion  Layer 
7:  Cathode  Gas  Channel 


Fig.  1.  Schematic  of  alkaline  fuel  cell. 


Table  1 

Transport  equations  (p\  liquid  phase,  y:  gas  phase,  cr:  solid  phase). 

Mass 

^PpPp  +  V -epPpVp=S™y 
ftSyfiy  +  V  ■  SypyVy  =  ~SJy 

Momentum 

^t{£pPpVp)  +  V  •  (spPpVp vp)  =  -SpVpp  +  V  •  (fif  Wvp)  +  epppg  +  SJ 

|( SyPyVy )  +  V  •  ( SyPyVyVy )  =  SyVpy  +  V  •  (/if  VVy)  +  SyPyg  +  SJ 

Species 

| ^pPpYpi)  +  V  ■  (epppvfy)  =  V  ■  (ppDf.  AYpi)  -^ip-  Vtp  +  S*.;  i  = 
H2,02,H20,  OH" 

UeyPyYyt)  +  V  •  (SyPyVyVy,)  =  V  •  (PyD^  AYy,)  +  i  =  H2 , 02  ,  H20 

Energy 

PCp‘:fr  f  +  (.CppepPpVp  +  CpyByPyVy )  ■  VT  = 

V  ■  (Xe&VT)  +  S^eact  +  Sphase  +  Sjoule 
Charge 

V  ■  (xfV0p)  +  V  •  (k^V  In  Cpe)  = 


2.3.  Ideal  solution 

The  Nernst-Planck  equation  is  used  to  calculate  the  flux  of 
charge  species  in  electrolyte.  For  ideal  solution,  the  flux  of  species 
k  is  obtained  from  Eq.  ( 1 1 )  [  7  ] : 


Cathode: 


Jk  =  -PDkVYk  -  pYkukZkFV0  +  f>Ykv 


(11) 


02  +  2H20  +  4e~  ->■  40H“  (2) 

The  anode  and  cathode  half  cell  voltage  is  calculated  by 

Fa  =  Ua  +  Via  +  0a  (3) 

EC  =  Uc  +  TJc  +  0C  (4) 

where  subscripts  a  and  c  denote  anode  and  cathode,  respectively. 
The  overall  cell  voltage  is  obtained  from 


Fee//  —  F c  —  Fa  —  {Uc  —  Ua)  (rjc  —  rja)  ( 0C  —  0a ) 


(5) 


where  F  is  the  potential  of  the  solid  phase,  0  is  the  potential  of 
the  liquid  phase,  and  r\  is  the  activation  overpotential.  The  Uc  -  Ua 
term  in  Eq.  (5)  is  the  theoretical  equilibrium  cell  voltage  calculated 
from  the  Gibbs  free  energy;  rjc  -  r)a  is  the  activation  overpotential, 
and  0C  -  0a  is  the  potential  drop  due  to  the  concentration  overpo¬ 
tential  and  IR  losses.  The  half-cell  open  potential  is  calculated  by 
Nernst  equation: 


Ua  =  U°a  - 


uc  =  au  - 


RuT 


. naF_ 


RUT 


naF 


in  [*£ 


P 


o2 


C° 

OH“ 


(6) 


(7) 


where  U°  is  the  theoretical  open  cell  potential  evaluated  at  stan¬ 
dard  concentration  C°  at  1  atm  and  given  temperature  T;  sk  is  the 
stoichiometric  coefficient  of  species  k;  na  is  the  number  of  electron 
transfer  at  anode.  The  temperature  effects  on  theoretical  open  cell 
potential  are  calculated  following  [14]: 


ut  =  U%8  +  {T  -298A5)[ 


dU° 

dT 


298 


=  -0.823  -  (T  -  298.15)  x  8.360  x  10 
U °  =  0.4011  -  (T  -  298.15)  x  1.6816  x  10 


-4 


-3 


(8) 

(9) 

(10) 


The  REIS  of  the  equation  are  the  Fickian  diffusion,  migration, 
and  convective  terms,  respectively.  For  concentrated  solutions,  the 


Table  2 

Source  terms  of  governing  equations  (Table  1 ). 


Mass  transfer  rate  at  the  liquid-gas  interface  of  catalyst  layer 


Anode  :  S™  =  -a*  ^  (  £r  ~  )  mh2o  +  apyDpMl 


hHo Ph2 -cHo 


Cathode:  S™  = -a,* 


^P y  ~  “ PY  8, 

Liquid-phase  viscous  drag 
Catalyst  layers  : 

GDL : 


(  Wr  -  itr  )  Mh2o  +  apyDp, 02  °2P°^ 


Mh2 
—  Mn. 


cv  _  SP»PVP 

cv  _  SyJlyVy 

— 


Ky 

Hydroxide  ions  reaction  rate  (anode  catalyst  layer  and  cathode  catalyst  layer) 
SP,0H-  =  -  f0H -)apjn,  [anode (n  =  a),  cathode (n  =  c)] 

Hydrogen  reaction  rate  and  dissolvation  rate  (anode  catalyst  layer) 

SAh2  =  - +  af>rDP. h2  Mh2 

Oxygen  reaction  rate  and  dissolvation  rate  (cathode  catalyst  layer) 

cm  MOn  .  ,  ^  H0oP09-C09  ,, 

$p,o2  ~  ~4rapcrh  +  a^yDp  Q2  Ti  Mo2 

Hydrogen  dissolvation  rate  (anode  catalyst  layer) 

^H2-^,H2HH2PrCH2MH2 

Oxygen  dissolvation  rate  (cathode  catalyst  layer) 

S™O2=-apyDp,02H-^Mo2 

Water  evaporation  rate  (anode  catalyst  layer  and  cathode  catalyst  layer) 

cm  „  °Y’ H20  (  ps  Ph20  \  j.. 

■^y,H20  8g  \RUT  RUT  j  MH20 

Charge  generation  rate  (anode  catalyst  layer  and  cathode  catalyst  layer) 

S^a  =  -apain,  anode  (n  =  a),  cathode  (n  =  c ) 

Heat  source  due  to  Joule  heating  (anode  catalyst  layer  and  cathode  catalyst 
layer;  separator) 

SU  = 

Heat  source  due  to  reversible  and  irreversible  reaction  heat  (anode  and 
cathode  catalyst  layer) 

Sh  =  ap„  (i„r,-&TAs) 

Heat  source  due  to  phase  change  (anode  and  cathode  catalyst  layer) 

=  — ^H2o^y>H2° 
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Table  3 

Constitutive  Equations. 


Local  reaction  current  (A  m-2 ) 


Effective  viscosity  (kg  irr 1  s-1 ) 

Effective  diffusion  coefficient  (m2  s-1 ) 
Effective  electrical  conductivity  (S  irr 1 ) 

Effective  diffusion  conductivity  (Am-1 ) 

Permeability  (m2) 

Capillary  pressure  (Pa) 

Current  flux  (A  irr 2 ) 

Effective  thermal  conductivity  (W  irr 1  K-1 ) 
Specific  heat  (J  kg-1  K_1 ) 

Electrical  potential  (V) 


fe)  te) 


m  ^ 
Kf  =  TiKP 


D,eff  _ 

Kfi  ~  ' 


,.eff  _  £K,, 

P'Y  —  Xy  ^ y 

Dff  =  ^Dy, 


s3dj 


kr, nw  =  Sr3,  kr,w  =  (1  -  Srf 


K  = 

16kK(l~ef 

gs^(1-417Sr  _  2.l20Sr2  +  1.263Sr3),  Sr  = 

(. K/e)l'z  1  ~sim 

ip  =  -  AC^’^V  In  C^e 

heff=SpXp  +  +  £aA.a 

Cpeff=  (SpPfiCpp  +  SyPyCpy  +  S0p0Cp0 )  /  CT 

0“=O,  0ca  =  Ecell 


generalized  Stefan-Maxwell  equation  is  used  [7,15]: 

i:©(Li)-(^)(v“iK‘fii  (,2) 

1  =  1 

i±k 

where  the  chemical  potential  of  species  k  is  given  by 

/ik  =  li°k  +  RuT  In  a*.  (13) 

In  Eq.  (13),  /x®  is  the  standard-state  chemical  potential,  and  ak 
is  the  activity  of  species  k.  The  current  density  is  calculated  by 
summing  the  fluxes  of  all  charge  species: 


•  \^zkFjk 
Mk 


(14) 


The  electrolyte  solution  in  AFCs  consists  of  charge  species  (e.g., 
I<+,  OH-)  and  solvent  (H20).  Following  Newman  [7],  the  mass  flux 


of  species  k  is  calculated  by 
Jk=  (  ~pDmVYk  +  pYkv 

where  tk  is  the  transference  number 


The  current  density  is  calculated  by 

N 

i  =  -kW0  -  y^zkFDkWCk 

k= l 

where  /c  is  the  KOH  solution  electrical  conductivity 

k = jykF*ukck 

k=  1 

and  Dm  is  the  KOH  diffusion  coefficient 


(15) 

(16) 

(17) 

(18) 

(19) 


Fig.  2.  Flowchart  of  CFD  calculation. 
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i  /  A.cnr2 


Fig.  3.  Comparison  of  CFD  calculation  with  experimental  data  [24]. 


where  subscripts  *+’  and  denote  the  positive  and  negative  ions, 
respectively. 

For  ideal  solution,  the  Nernst-Einstein  equation  is  used  to  cal¬ 
culate  the  mass  diffusivity  of  species  k 


Dk  —  ukRuT  (20) 

where 


=  tkKe 
Z2kVkF2Ce 


2.4.  Non-ideal  solution 

For  non-ideal  solution,  the  mass  flux  of  species  k  is  calculated 
by 

A  =  ~  pDmVYk  +  pYkv  (21 ) 


CKoh  /M 

Fig.  4.  KOH  electrical  conductivity  as  a  function  of  KOH  molarity;  T  =  20,  60,  80, 
100°C. 


Table  4 

Limiting  ionic  equivalent  conductance  [16]. 


Temperature  (°C) 

I<+  (cm2  £2-1  mol-1 ) 

OH- 

(cm2  mol-1) 

25 

73.5 

198.3 

45 

103.4 

55 

119.2 

75 

152.9 

100 

195 

450 

125 

240 

where 


MD7°)K;tc) 

(22) 

D  D0+D0_(z+  -Z-) 

(23) 

z+D0+  —  z_D0_ 

CT  =  C+  +  C_ 

(24) 

pot?  +  P- 
+  =  P 

(25) 

In  Eqs.  (22)-(25),  D0+  and  D0_  are  the  binary  diffusion  coef¬ 
ficient,  Ce  is  the  KOH  concentration,  CT  is  the  total  molar 
concentration,  y±iC  is  the  mean  molar  activity  coefficient,  and  is 
the  transference  number.  The  dependency  of  transference  number 
on  electrolyte  concentration  is  calculated  by 


t°  X°+  + 


(26) 


where  X9  is  the  limiting  ionic  equivalent  conductance  of  ion  i.  The 
values  of  A.9  at  several  selected  temperatures  are  given  in  Table  4 
[16].  The  current  density  is  calculated  from 


i  =  -/cV0  -  /cDV  In  Ce 


(27) 


where  kd  is  the  diffusion  conductivity: 


d  lny±;C\ 
d  In  Ce  j 


(28) 


The  non-ideal  solution  effects  are  calculated  from  the  non-unity 
activity  coefficients  that  appear  in  the  species  diffusion  coefficients 
(Eq.  (22))  and  diffusion  conductivity  (Eq.  (28)). 

The  primary  and  secondary  (shunt)  currents  are  calculated  by 
the  charge  conservation  equation.  The  source  term  in  the  charge 
equation  is  used  to  calculate  the  charge  transfer  between  the  solid 
and  liquid  phase.  The  Butler-Volmer  equation  is  used  to  calculate 
the  H2  oxidation  reaction  (HOR)  and  02  reduction  reaction  (ORR) 
rates.  An  infinitely  large  electrical  conductivity  is  assumed  for  the 
solid  phase.  Namely,  the  electrode  is  assumed  to  have  an  equipo- 
tential  surface.  Conversely,  @s  =  0  and  @s  =  Eceu  are  set  for  the  anode 
and  cathode,  respectively. 

A  direct  two-electron  transfer  reaction  is  assumed  for  HOR 
and  a  four-electron  transfer  reaction  for  ORR,  e.g.,  Eqs.  (1)  and 
(2),  respectively.  The  exchange  current  density  is  assumed  to  be 
5.0  x  1 0-4  A  cm-2  for  HOR  [17],  and  5.0  x  10-8  A  cm-2  for  ORR  [18]. 
The  specific  catalyst-electrolyte  interface  area,  is  estimated  fol¬ 
lowing  Jo  and  Yi  [13].  The  specific  gas-electrolyte  interface  area, 
ag,  and  the  thickness  of  the  electrolyte  film,  8i,  are  specified  using 
the  reported  values  1 05  m_1  and  1 0_1  p,m,  respectively  [  1 9,20].  The 
thickness  of  the  GDL  characteristic  length,  8g,  is  set  to  10  pm  [21  ]. 
The  entropy  changes  of  HOR  and  ORR  are  taken  from  [22].  The 
electrochemical  parameters  used  in  the  calculation  of  the  baseline 
condition  are  summarized  in  Table  5.  The  geometric  parameters, 
along  with  the  boundary  conditions,  are  given  in  Table  6.  Thermo¬ 
dynamic  properties  used  in  the  calculation  are  listed  in  Table  7. 
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Fig.  5.  KOH  diffusion  coefficient  as  a  function  of  KOH  molarity;  T=20,  60,  80, 100  °C;  (a)  non-ideal  solution  formulation,  and  (b)  ideal  solution  formulation. 


3.  Solution  method 

3.1.  Computational  domain 

The  computational  domain  is  divided  into  seven  regions 
(domains)  as  illustrated  in  Fig.  1.  Domains  1,2,  and  3  are  the  anode 
gas  channel,  gas  diffusion  layer  (GDL),  and  catalyst  layer,  respec¬ 
tively.  Domain  4  is  the  separator.  Domains  5,  6,  and  7  are  the 
cathode  catalyst  layer,  GDL,  and  gas  channel,  respectively.  H2  is 
transported  from  Domain  1  to  Domain  2  by  convection  and  diffu¬ 
sion,  and  dissolves  in  KOH  solution.  The  dissolved  H2  reacts  with 
hydroxide  ions  (OH-)  and  forms  H20  at  the  anode  catalyst  layer 
(Domain  3).  Electrons  are  released.  The  electrochemical  reaction  at 
the  cathode  catalyst  layer  (Domain  5)  involves  returned  electrons 
from  external  circuit,  dissolved  02,  and  H20. 02  is  transported  from 
Domain  7  to  Domain  6,  and  dissolves  in  KOH  solution.  The  OH-  ions 
are  produced  in  Domain  5,  carried  away  by  diffusion,  convection 
and  migration  actions,  and  consumed  in  Domain  3.  The  electrolyte 


flows  through  the  separator  (Domain  4);  it  also  takes  away  the  reac¬ 
tion  heat  and  H20  produced  in  Domain  3.  Two  extended  channels 
(region:  DEFG  and  OPQN)  at  separator  inlet  (edge:  DG)  and  exit 
(edge:  QN)  are  implemented  to  capture  the  external  electrolyte 
flow. 

3.2.  CFD 

FLUENT®  was  used  to  solve  the  coupled  governing  equations. 
The  flow  chart  of  the  calculation  procedure  is  illustrated  in  Fig.  2. 
The  user  defined  functions  (UDFs)  were  written  to  calculate  the 
physical  properties  and  the  source  terms  for  solving  the  governing 
equations,  as  well  as  to  satisfy  the  prescribed  boundary  conditions. 
The  graphical  user  interface  (GUI)  was  used  to  input  the  UDFs  to  cor¬ 
responding  governing  equations  and  the  user  defined  scalar  (UDS) 
was  introduced  to  solve  the  charge  equation.  Two  convergence  cri¬ 
teria  were  enforced  to  ensure  that  (a)  the  local  cell  voltage  is  the 
same  everywhere  along  the  separator,  and  (b)  the  total  current  den- 


Fig.  6.  Reactant  solubility  in  KOH  solution  as  a  function  of  KOH  molarity,  T  =  20,  60,  80, 100  °C,  p  =  1  atm;  (a)  H2  and  (b)  02. 
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Table  5 

Baseline  conditions:  electrochemical  kinetic  parameters  for  the  anode  and  cathode. 
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Anode 

Cathode 

Parameter 

Value  (unit) 

Parameter 

Value  (unit) 

na 

2 

nc 

4 

Sh2 

1 

So2 

-1 

SOH- 

2 

SOH- 

4 

Sh20 

-2 

Sh20 

-2 

Qh2 

1 

Qo2 

1 

Qoh~ 

2 

Qoh- 

4 

<Zh2o 

2 

Qh2o 

2 

vK+ 

1 

vK+ 

1 

VOH~ 

1 

VOH~ 

1 

zK+ 

1 

zK+ 

1 

ZOH~ 

-1 

ZOH“ 

-1 

OLaana 

1.5 

acanc 

2.5 

aacna 

0.5 

accnc 

1.5 

*o 

5.0  x  10“4  (A cm-2)3 

'o 

5.0  xlO-8  (Acrrr2)b 

af 

1.0  x  107  (m-1)' 

of 

1.0  xlO7  (m-1)0 

al 

6.0  x  105  (irr1)d 

ag 

6.0  x  105  (m_1)d 

1.0  x  10“7  (m)e 

sc, 

1.0  xl0“7  (m)e 

n 

1.0  x  10“5  (m)f 

% 

1.0  x  10“5  (m)f 

C»2 

6.0913  x  10“7  (mol cm-3) 

Co2 

3.4905  x  10“7  (mol cm-3) 

ASa 

161.2(J  moF1  K-1  )s 

A  Sc 

-648.0  (J  moF1  K-1  )g 

a  Tilaketal.  [17]. 
b  Kinoshita  [18]. 
c  Jo  and  Yi  [13]. 
d  Kenjo  [19]. 
e  Li  et  al.  [20]. 
f  Nam  and  Kaviany  [21  ]. 
g  Lampinen  and  Fomino  [22]. 


Table  6 

Baseline  conditions:  AFC  structural  parameters  (I:  thickness;  H:  height)  and  bound¬ 
ary  condition. 


Parameter 

Value  (unit) 

Parameter 

Value  (unit) 

t separator 

300  ((Jim) 

£ separator 

1.0 

La 

GDI 

250  (|Jim) 

Sa 

CGDL 

0.7 

Lc 

GDI 

250  (|Jim) 

£c 

CGDL 

0.7 

Kat 

25  (|Jim) 

£a 

bcat 

0.7 

LCcat 

25  (|Jim) 

£c 

bcat 

0.7 

La 

channel 

1  (mm) 

tgdl 

1.2 

L channel 

1  (mm) 

tc 

lGDL 

1.2 

H separator 

50  (mm) 

ra 

Lcat 

1.2 

jjextended 

1 1  separator 

50  (mm) 

TC 

1 cat 

1.2 

Inlet  condition:  Ce  =  7  M,  T=  80  °C,  p  =  4.1  atm;  v a,  vc,  ve  =  0.2,  0.1 ,  0.01  m  s”1 ,  RH  =  0%. 
Wall:  adiabatic  and  impermeable,  9</>/9n  =  0,  ((f)  =  T,  Yi,  0). 

Electrolyte  inlet  and  exit:  9d>/9n  =  0. 


Fig.  7.  Comparison  of  non-ideal  solution  versus  ideal  solution  results:  polarization 
curve,  isothermal  boundary  condition. 


sity  equals  that  calculated  by  the  specified,  averaged  cell  current 
density. 

4.  Results  and  discussion 

4.1.  Validation 

The  CFD  results  of  baseline  condition  (i.e.,  parameters  spec¬ 
ified  in  Tables  5  and  6)  were  validated  by  experimental  data 
reported  in  [24].  The  catalysts  used  in  [24]  were  PtPd  (loading  at 
1 0  mg  cm-2 )  and  AuPt  (loading  at  20  mg  cm-2 )  for  anode  and  cath¬ 
ode,  respectively.  The  computation  was  performed  to  represent  the 
operating  condition:  pure  02  and  pure  H2,  80 °C  and  4.1  atm,  and 
KOH  concentration  of  7M.  Excellent  agreement  (discrepancy  less 
than  5%)  was  seen  for  current  density  greater  than  0.30  A  cm-2, 
and  good  agreement  (discrepancy  less  than  10%)  for  current  den¬ 
sity  below  0.30  A  cm-2 .  The  comparison  is  shown  in  Fig.  3.  It  should 
be  noted  that  estimates  were  made  on  the  parameters  that  were 
not  available  in  the  literature.  These  parameters  are  the  effective 
liquid-solid  surface  area,  the  effective  liquid-gas  surface  area,  and 
the  liquid  diffusion  film  thickness.  The  calculation  also  showed  that 
the  cathode  overpotential  was  the  predominant  loss  mechanism 
that  accounted  for  about  70%  of  the  potential  losses  [27]. 

4.2.  Property  effects 

The  effects  of  thermodynamic  and  transport  properties  on  AFC 
operation  were  examined.  The  electrolyte  conductivity  directly 
impacts  the  IR  losses.  The  reactant  dissolving  rates  can  set  the  lim¬ 
iting  current  density  condition.  The  KOFI  conductivity  varies  with 
temperature  and  molarity,  e.g.,  see  Fig.  4.  At  80  °C,  the  conductiv¬ 
ity  peaks  at  around  7  M,  which  is  the  molarity  condition  chosen 
for  AFC  when  operated  at  80  °C  [24].  Deviation  from  this  molarity, 
for  example,  due  to  the  presence  of  concentration  gradient,  will 
increase  the  IR  losses  and  decrease  the  cell  voltage. 

The  reactant  (i.e.,  H2  or  02 )  dissolving  rates  are  functions  of  reac¬ 
tant  solubility,  diffusion  coefficients,  specific  gas-liquid  interfacial 
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nonideal  solution 
T=  80*C,  i=0.8  A/cm1 


ideal  solution 
T=  80°C,  i=0.8  A/cm2 


X/m 


X/m 


Fig.  8.  Comparison  of  non-ideal  solution  versus  ideal  solution  results:  isopleths  of  KOH  concentration  in  separator  and  extended  channel;  i  =  0.8  Acm-2;  anode  and  cathode 
located  at  x  =  0  and  0.001 6  m. 


6  6 


Fig.  9.  Streamwise  profiles  at  cathode-electrolyte  interface  of  non-ideal  solution  formulation,  i  =  1 .5  A  cm-2 ;  <5  is  the  dimensionless  distance,  inlet  at  8  =  0  and  exit  at  d  =  1  and 
exit;  (a)  temperature,  (b)  KOH  molarity,  (c)  oxygen  partial  pressure,  and  (d)  dissolved  oxygen  molarity. 
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nonideal  solution 
adiabatic,  i=1.5  A/ cm2 

0.15 


0.10 


a 


0.05 


0.0014  0.0016 


nonideal  solution 
T=  80®C,  i=1.5  A/cm2 

0.15n - 


c^oh  m 


■ 


7.93 

7.79 

7.66 

7.52 

7.38 

7.24 

7.10 

6.97 

6.83 

6.69 

6.55 

6.41 

6.28 

6.14 

6.00 


0.00 


j _ . _ I _ . 

0.0014 


0.0016 


X/in  X/rn 

Fig.  10.  Effects  of  adiabatic  versus  isothermal  boundary  condition  on  calculated  isopleths  of  KOH  concentration  in  separator  and  extended  channel;  i=  1.5  A  cm-2,  non-ideal 
solution  formulation;  anode  and  cathode  located  atx  =  0  and  0.0016  m. 


areas,  and  liquid  film  thickness.  At  constant  temperature,  the  KOH 
diffusion  coefficient  reaches  its  minimum  over  the  concentration 
range  0-1  M,  e.g.,  see  Fig.  5(a).  The  diffusion  coefficient  increases 
with  increasing  temperature  and  concentration.  It  increases  by  a 
factor  of  6  when  the  temperature  is  increased  from  20  to  1 00  °C,  but 
the  increase  is  less  than  20%  when  the  concentration  is  increased 
from  0  to  12  M.  Also  plotted  in  Fig.  5  is  the  calculated  species  dif¬ 
fusion  coefficient  when  the  ideal  solution  assumption  is  invoked, 
i.e.,  Fig.  5(b)  for  temperatures  set  to  20,  60,  80,  or  1 00  °C.  The  com¬ 
parison  shows  that  when  the  ideal  solution  assumption  is  invoked 
it  under-estimates  the  KOH  diffusion  coefficient  and  predicts  an 
opposite  trend  with  increasing  KOH  concentration  when  compared 
to  the  non-ideal  solution  results. 

Solubility  of  H2  and  02  in  KOH  solution  is  of  the  same  order 
of  magnitude  as  shown  in  Fig.  6(a)  and  (b),  respectively.  Solubility 
decreases  by  about  2  orders  of  magnitude  with  increasing  molarity 
over  the  range  0-12  M,  but  the  change  is  much  less  with  increas¬ 
ing  temperature  over  the  range  20-1 00  °C.  In  AFC,  KOH  is  produced 
at  cathode  and  consumed  at  anode,  resulting  in  higher  concentra¬ 
tions  at  cathode  and  lower  at  anode.  Reaction  therefore  favors  H2 
dissolution  into  KOH  solution,  and  retards  02  dissolution. 

The  presence  of  KOH  concentration  gradient  across  the  elec¬ 
trolyte  increases  the  IR  losses  when  the  concentration  deviates 
from  the  design  point  of  7  M.  To  reduce  the  concentration  gradient, 
one  can  increase  the  rate  of  KOH  transport  from  cathode  to  anode, 
or  increase  the  KOH  flow  rate.  For  a  fixed  fuel  cell  configuration,  the 
former  can  be  achieved  by  operating  the  fuel  cell  at  higher  tem¬ 
peratures  to  increase  the  diffusion  coefficient  (e.g.,  see  Fig.  5(a)). 
However,  increasing  the  temperature  will  also  increase  the  H20 


evaporation  rate,  which  in  term  will  increase  the  KOH  concentra¬ 
tion  and  reduce  02  solubility  at  cathode.  The  latter  can  be  achieved 
at  the  expense  of  parasitic  losses  associated  with  the  increase  of 
the  pumping  power  to  re-circulate  the  electrolyte. 

4.3.  Ideal  solution  assumption 

Isothermal  boundary  condition  was  chosen  to  examine  the 
effects  of  ideal  solution  assumption  on  the  calculation  of  polariza¬ 
tion  curves.  Fig.  7  shows  that  when  the  ideal  solution  assumption 
is  relaxed  the  calculated  limiting  current  density  is  extended  from 
1.0  to  1.8  A  cm-2,  which  is  a  result  of  the  increased  diffusive  trans¬ 
port  of  KOH  in  the  electrolyte.  As  shown  by  the  isopleths  illustrated 
in  Fig.  8,  the  ideal  solution  assumption  results  in  concentrated 
KOH  at  cathode  (in  the  range  of  5.3-9.5  M  for  calculation  with  i 
set  to  0.8  A  cm-2)  and  reduced  KOH  concentration  at  anode  (in 
the  range  of  5.5-7.0M).  The  non-ideal  solution  calculation,  how¬ 
ever,  yields  a  more  uniform  KOH  concentration  (in  the  range  of 
6.4-7. 7  M).  These  results  can  be  explained  by  comparing  the  cal¬ 
culated  diffusion  coefficients  given  in  Fig.  5(a)  and  (b).  The  value 
of  diffusion  coefficients  calculated  based  on  the  ideal  solution  for¬ 
mulation,  Fig.  5(b),  was  lower  than  that  based  on  the  non-deal 
solution  formulation,  Fig.  5(a).  The  diffusive  transport  of  KOH  from 
cathode  to  anode  was  lower  with  the  ideal  solution  formulation. 
Consequently,  higher  KOH  concentrations  were  predicted  at  cath¬ 
ode,  which  in  term  reduced  the  02  dissolving  rate  and  lowered  the 
limiting  current  density.  As  noted  earlier,  the  maximum  electrical 
conductivity  occurs  at  KOH  concentration  around  7  M.  A  uniform 
KOH  concentration  around  7  M  is  expected  to  result  in  higher  cell 
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Table  7 

Thermodynamic  properties. 

Electrical  conductivity  (Tin  K;  range:  up  to  9  M  and  above  0°C)  [16]: 

Ke  (S  cm-1 )  =  0.0262 W+  6.7  x  10-4  W(T-  273.15)  -  4.8  x  10 -4W2  -  8.8  x  10-6  W2  ( T - 273.15) 

Diffusion  coefficient  (Ce  in  mol  cm-3,  Tin  K;  range:  up  to  14  M  and  90 °C)  [28]: 

Dm  (cm2  s”1 )  =  exp(-16.489  +  0.02015T  -  8.1607Q0  5  +  286.2Ce  -  2539.8C1-5  +  7207.5 C2) 

Density  (p  in  kgm-3,  m  in  molkg-1  and  Tin  K;  range:  up  to  14M  and  100°C)  [29,30]: 
p  =  (2.1 1942  x  103  +  1.03561  x  102  m- 3.67252  m2  +  9.30540  x  10-2  m3  -  1.11476  x  10-3  m4)  x  T-013255 
Activity  coefficient  (Tin  K  range:  m  from  2  to  18  molkg-1  and  T from  0  to  200 °C)  [31]: 
y±tC  =  Y±,m(Pu2o/(pe  -  MeCe));  log  y±,m  =  473.520  -  76.168/T  -  8.1727T  +  0.05818T2  -  2.1864  x  10“4T3  +  4.5717  x  10“7T4  -  5.0464  x  10“10T5  +  2.2960  x 
1013T6  - (0.079599 -51.367/T)  m  +  (0.44109  -  88. 190/T)  logm 
Specific  heat  (T in  K;  range:  25-150°C)  [32]: 

Cp  (J  kg-1  K-1 )  =  2358.948  +  2.0934T 

Thermal  conductivity  (Ce  in  molL-1;  a\  =-7.19  x  10-4,  02  =  -1.0  x  10-4;  range:  1-16 M)  [33]: 

A.38°c  (kcal  m“1h“1k_1)  =  0.538 +  aiCe  +  a2Ce2;  (AA38°c)soin  =  (A A38°c)h2o 
A.h2o  (W m_1  K_1 )  =?  1.487188  x  10“8T3  -  1.143515  x  10“5T2  +  2.193975  x  10“3T  +  0.559401 ;  Tin°C 
Dynamic  viscosity  (Tin °C;  range:  3-11  M  and  10-60°C)  [34]: 

IL( kgm-1  s-1)  =  (1.71468  x  10-2  -  1.02529  x  10-1  W+  2.72493  x  10-1  W2)T-0-5414 
Gas  solubility  in  electrolytes  ( kscx  in  m3  kmol-1;  range:  25-100 °C)  [12,13,35]: 

log (xJVx,-)  =  ksacCe]  kscx, h2  =  0.129;  kscx, 02  =  -4.1315  x  10“8T3  +  1.4405  x  10“5T2  -  1.5914  x  10“3T  + 2.1139  x  10”1 

Ck  =  Hekpk\  Hek  =  (C0  +  2Ce)(xk/(l  -  xfc));  C0  :  solvent,  Ce  :  salt 
Gas  solubility  in  pure  H20  ( H  in  kmol  m-3  atm-1 ;  Henry’s  law)  [36]: 
logH*  =  -  (1.142  -  2.846(1 /T)*  +  2.486(l/T)*2  -  0.9761  (1/T)*3  +  0.2001  (1/T)*4) 

H *  =  H/Hmax  ;  ( 1  IT)*  =  ( 1/T—  1  /Tr )/( 1  /Tmax  -  1  /Tr) ;  Tr  =  647  K  (critical  temperature  of  water) 

Hmax  x  1 0-4  =  7.08  (02 ),  7.54  (H2),  12.39  (N2) ;  (1/Tmax)  x  103  =  2.73(02),  3.09(H2),  2.80(N2) 

Diffusion  coefficient  of  dissolved  gases  (D  in  cm2  s-1 ;  range:  KOH  up  to  15  M,  25-100  °C)  [37]: 

DlH  =  (1.3004  x  10~6  -  8.4347  x  10“6Y°-5  +  2.7822  x  10“5Y- 4.1593  x  10“5Y15  +  2.2415  x  10“5y2)T1  28968 
DlQ2  =  (5.5421  x  10~7  -  2.3393  x  10“6Y°-5  +  6.3923  x  10“6Y  -  9.6313  x  10“6Y15  +  5.4576  x  lO^Y^T1-20838 
Vapor  pressure  (p  in  bar,  T in  K;  range:  0-300 °C  KOH  up  to  25  molkg-1)  [38,39]: 

log  pw  =  a  +  b  logp6  ;  a  =  (-1.508m  -  1.6788m2  +  2.25887  x  10“3m3)  x  10“2,  b  =  1  -  (1.2062m  +  .56024m2  -  7.8228  x  10“3m3)  x  10“3 
log  p6  =  35.4462  -  3343.93 /T  -  10.91  log  T  +  .0041645T 
Electrode  thermal  conductivity  [25,26]: 

Xa  =  Zc  =  1.5  Wm-1  K1 


voltages  when  operated  in  the  IR  dominated  regime,  e.g.,  see  Fig.  8 
for  the  results  of  current  density  in  the  range  0.4-0.8  A  cm-2.  The 
property  considerations  suggest  that  maintaining  a  uniform  KOH 
concentration  in  the  electrolyte,  for  example  at  7  M,  be  a  design 
strategy  to  extend  the  limiting  current  density  condition,  and  to 
reduce  the  IR  losses  in  AFC. 

4.4.  Thermal  boundary  conditions 

To  examine  the  effects  of  thermal  conditions  on  the  calculation 
of  polarization  curves,  the  baseline  condition  of  isothermal  bound¬ 
ary  condition  (at  80  °C)  was  changed  to  one  subject  to  adiabatic 
boundary  condition.  One  major  effect  of  the  adiabatic  boundary 
condition  was  that  it  increased  the  cathode  temperature  along  the 
streamwise  direction  from  82  °C  at  inlet  (5  =  0)  to  just  above  1 00  °C 
at  exit  (5  =  1),  see  Fig.  9(a).  As  discussed  earlier,  02  solubility  is  rel¬ 
atively  insensitive  to  temperature  change  over  the  range  studied, 
cf.,  Fig.  6(b);  however,  the  dissolved  02  level  was  orders  of  magni¬ 
tude  higher  at  the  exit,  Fig.  9(d).  The  oxygen  partial  pressure  was 
about  two  orders  of  magnitude  lower  at  the  exit,  e.g.,  see  Fig.  9(c). 
The  higher  dissolved  02  level  was  due  to  the  higher  KOH  diffusion 
coefficient  calculated  by  the  adiabatic  boundary  condition.  The  dif¬ 
fusion  coefficient  increased  by  an  order  of  magnitude  when  the 
temperature  was  increased  from  80  to  100  °C,  cf.  Fig.  5(a).  The  KOH 
transport  rate  from  cathode  to  anode  was  increased  accordingly, 
and  it  resulted  in  lower  KOH  concentrations  at  the  cathode,  Fig.  9(b). 
The  isopleths  of  the  results  subject  to  adiabatic  and  isothermal 
boundary  conditions  are  illustrated  in  Fig.  10.  For  isothermal  con¬ 
dition,  the  KOH  concentration  remained  high  along  the  cathode 
(i.e.,  along  the  y-direction  at  x  =  0.001 6  m)  and  low  along  the  anode 
(i.e.,  along  the  y-direction  at  x  =  0.00  m).  The  growth  of  the  concen¬ 


tration  layer  was  reminiscent  to  that  of  the  boundary-layer  flow. 
However,  this  was  not  the  case  for  the  adiabatic  boundary  condi¬ 
tion.  The  increase  of  diffusion  coefficient  along  the  cathode  helped 
to  transport  away  KOH  to  anode.  The  thickness  of  the  concentra¬ 
tion  boundary  decreased,  the  KOH  concentration  decreased,  and 
dissolved  02  level  increased.  The  net  effects  are  that  the  adiabatic 
boundary  condition  extended  the  limiting  current  density  from 
1.8  to  2.3  A  cm-2  (cf.,  Fig.  11).  It  is  also  noted  that  both  thermal 


Fig.  11.  Effects  of  adiabatic  versus  isothermal  boundary  condition  on  calculated 
polarization  curve;  non-ideal  solution  formulation. 


G.  Zhou  et  al.  /  Journal  of  Power  Sources  196(2011)  4923-4933 


4933 


conditions  yielded  similar  cell  voltages  for  i<  1.6  A  cm-2,  although 
the  adiabatic  boundary  condition  predicted  a  significantly  higher 
limiting  current  density. 

5.  Summary  and  conclusion 

A  computational  investigation  was  conducted  that  relaxed 
the  ideal  solution  assumption,  and  accounted  for  thermodynamic 
solubility  of  the  reactants  and  property  variations  due  to  the  tem¬ 
perature  effects.  For  the  condition  studied,  at  which  the  cathode 
overpotential  accounts  for  about  70%  of  the  cell  potential  losses, 
the  following  conclusions  can  be  drawn: 

1.  The  ideal  solution  formulation  underestimates  the  limiting  cur¬ 
rent  by  about  50%.  This  level  of  discrepancy  suggests  that  the 
ideal  solution  assumption  should  not  be  applied  in  the  modeling 
and  simulation  of  AFC  for  conditions  similar  to  that  considered 
in  this  study. 

2.  The  thermal  condition  can  significantly  change  the  results  of  the 
calculation.  The  adiabatic  boundary  condition  can  increase  the 
limiting  current  by  30%.  Conversely,  thermal  management  can 
significantly  extend  the  limiting  current  condition,  thus  increase 
the  power  density  of  the  fuel  cell. 

3.  The  computational  results  showed  that  a  nearly  uniform  distri¬ 
bution  of  KOH  in  the  electrolyte  was  responsible  for  the  higher 
limiting  current  density  condition  calculated  in  the  study.  Con¬ 
versely,  a  design  strategy  to  yield  a  uniform  KOH  concentration 
in  the  electrolyte  (e.g.,  at  the  design  point  of  7  M)  will  be  an  effec¬ 
tive  measure  to  increase  the  limiting  current  density  condition. 
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